Validation round ensembled model scoring

Introduction

Often times, an ensemble of methods will perform better than the individual method. This known as the “wisdom of the crowds” phenomenon. An easy way to generate an ensemble prediction is to take the mean, median, or weighted average of all of the predictions. You then can score this “prediction” as you would any other prediction file to assess it’s performance relative to the submissions.

Another consideration is that the wisdom of the crowds method sometimes applies only to a certain point. That is, if you order all of the submitted predictions from high to low performance, there may be a point after which you no longer want to add a prediction to your ensemble method. A good visualization of this can be found in Supplemental Figure 8 here, where you can see how the ensemble score changes after adding additional predictions. Performance peaks with an ensemble of the top four predictions, but does not improve with additions of further models.

A final consideration here is that certain models may be better at predicting certain components of the problem. You may be able to more strategically ensemble methods if you assess the performance of the submitted models on categorical subsets of the data (e.g. for a method predicting among four cancer subtypes, some methods may be better at predicting subtype 1 vs the other methods, while other methods may be better at predicting subtype 3 than other methods). Weighting the predictions from these methods strategically may result in better ensemble performance.

First, import packages, scoring functions, and challenge data.

library(tidyverse)
library(reticulate)
# library(reactable)
use_condaenv("synapse-2") #conda environment with synapse >2.0 installed
synapse <- import('synapseclient')
syn <- synapse$Synapse()

Scoring helper functions.

# 
calculate_weight <- function(x){
  weight <- dplyr::if_else(x == 0, 1,
            dplyr::if_else(x == 1, 2,
            dplyr::if_else(x>=2 & x<=3, 2.143547,
            dplyr::if_else(x>=4 & x<=7, 3.863745,
            dplyr::if_else(x>=8 & x<=20, 8,
            dplyr::if_else(x>=21 & x<=55, 16,
            dplyr::if_else(x>=56 & x<=148, 32,
            dplyr::if_else(x>=148, 64, 0))))))))
  return(weight)
}

# 
sum_weighted_error <- function(gold, pred, weight){
  sum(weight*abs(gold - pred))
}

# 
rmse <- function(gold, pred){
  sqrt(mean((gold - pred) ** 2))
}

# Scoring function for joint weighted sum RMSE.
score_weighted_rmse <- function(df, pred) {
  df$pred_score <- pred
  df %>% group_by(Patient_ID, weight) %>%
    summarize(patient_rmse = rmse(log2(score+1), log2(pred_score+1))) %>% 
    ungroup() %>% 
    summarize(result= sum(weight*patient_rmse)/sum(weight)) %>%
    pluck("result")
}

Goldstandard

For the validation round, we have two expert reader measurements. Previous analysis showed that the validation predictions scored relatively concordantly between these two gold-standards, so to simplify this analysis we’ll take the average of the two datasets for each measurement.

The assessment of subchallenge 1 utilizies known SvH scores and subchallenges 2 & 3 known individual joint narrowing and erosion scores. After adding a weight for each joint score, split the goldstandard accordingly.

gold.sc1 <- gold %>% 
  select(Patient_ID, Overall_Tol, weight)

Error: Can't subset columns that don't exist.
x Column `weight` doesn't exist.
Run `rlang::last_error()` to see where the error occurred.

---

Bootstrap Submissions

Read in prediction files, combine, then bootstrap the predictions + a gold standard 1000 times to calculate 1000 scores per prediction.

p

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

p

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

bs_indices.sc2 <- matrix(1:nrow(ensemble.sc2), nrow(ensemble.sc2), N) %>%
  apply(2, sample, replace = T)
boot.sc2 <- apply(bs_indices.sc2, 2, function(ind) {
  tmp.gold <- ensemble.sc2[ind, c(1:4)]
  apply(ensemble.sc2[ind, -c(1:4)], 2, function(pred) {
    score_weighted_rmse(tmp.gold, pred)
  })
}) %>%
  t()
p

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

bs_indices.sc2 <- matrix(1:nrow(ensemble.sc2), nrow(ensemble.sc2), N) %>%
  apply(2, sample, replace = T)
boot.sc2 <- apply(bs_indices.sc2, 2, function(ind) {
  tmp.gold <- ensemble.sc2[ind, c(1:4)]
  apply(ensemble.sc2[ind, -c(1:4)], 2, function(pred) {
    score_weighted_rmse(tmp.gold, pred)
  })
}) %>%
  t()
### SC3
submissions.sc3 <- lapply(names(pred_filenames), function(team) {
  read.csv(pred_filenames[[team]]) %>% 
    select(-Overall_Tol) %>%
    gather('joint', !!team, -Patient_ID)
}) %>%
  reduce(left_join, by=c("Patient_ID", "joint")) %>%
  left_join(gold_joints, by=c("Patient_ID", "joint"))
results.sc3 <- submissions.sc3[, c(
    "Patient_ID", "joint", "weight", "score", (query %>% arrange(sc3))$submitterid
  )] %>%
  filter(grepl(".+_E__.+", joint))
ensemble.sc3 <- sapply(5:ncol(results.sc3), function(x){
  new_col <- results.sc3 %>% 
    select(5:x) %>% 
    mutate(!!sym(as.character(x)) := rowMeans(across(where(is.numeric)))) %>% 
    pluck(as.character(x))
}) %>% 
  as.data.frame()
ens_names <- c(glue::glue("{colnames(results.sc3)[5]}"), glue::glue("+{colnames(results.sc3)[-1:-5]}"))
colnames(ensemble.sc3) <- ens_names
ensemble.sc3 <- ensemble.sc3 %>% 
  add_column(weight = results.sc3$weight,
             score = results.sc3$score,
             Patient_ID = results.sc3$Patient_ID,
             joint = results.sc3$joint)
ensemble.sc3 <- ensemble.sc3[,c("Patient_ID", "joint", "weight", "score", ens_names)]
p<-lapply(ens_names, function(x){
  cor <- cor(ensemble.sc3$score, ensemble.sc3[[x]], method = 'spearman') 
  ggplot(data = ensemble.sc3) +
    geom_point(aes(x = score, y = !!sym(x))) +
    ggtitle(glue::glue("spearman = {cor}"))
})
p

[[1]]

[[2]]

[[3]]

[[4]]

[[5]]

[[6]]

[[7]]

[[8]]

[[9]]

[[10]]

[[11]]

[[12]]

[[13]]

[[14]]

bs_indices.sc3 <- matrix(1:nrow(ensemble.sc3), nrow(ensemble.sc3), N) %>%
  apply(2, sample, replace = T)
boot.sc3 <- apply(bs_indices.sc3, 2, function(ind) {
  tmp.gold <- ensemble.sc3[ind, c(1:4)]
  apply(ensemble.sc3[ind, -c(1:4)], 2, function(pred) {
    score_weighted_rmse(tmp.gold, pred)
  })
}) %>%
  t()

---

Bayes Factor Analysis

Use our challengescoring package to compute Bayes factors using a matrix of scores, setting the refPredIndex as the number of the column that contains the top prediction (the reference prediction).

In addition to computing BF against the top performer, we are also interested in comparing against all columns. Some tweaking to the computation will be required, as the reference index may no longer be the best, which in consequence, means the inversion of K will be dependent on which submission it is being compared with.

library(challengescoring)
computeBayesFactorWhereRefIsNotBest <- function(bootstrapMetricMatrix,
                               refPredIndex,
                               invertBayes){

    M <- as.data.frame(bootstrapMetricMatrix - bootstrapMetricMatrix[,refPredIndex])
    K <- apply(M ,2, function(x) {
      k <- sum(x >= 0)/sum(x < 0)
      if(sum(x >= 0) > sum(x < 0)){
      return(k)
      }else{
      return(1/k)
      }
    })
    if(invertBayes == T){K <- 1/K}
    K[refPredIndex] <- 0

    return(K)
}

# Top performer: Team Shirin
bayes_top.sc1 <- computeBayesFactorWhereRefIsNotBest(boot.sc1, refPredIndex=1, invertBayes=F) %>%
  as_tibble(rownames = "submission") %>%
  rename(bayes = value)

# Top performer: Hongyang Li and Yuanfang Guan (column 1)
bayes_top.sc2 <- computeBayesFactorWhereRefIsNotBest(boot.sc2, refPredIndex=1, invertBayes=F) %>%
  as_tibble(rownames = "submission") %>%
  rename(bayes = value)

# Top performer: Gold Therapy (column 1)
bayes_top.sc3 <- computeBayesFactorWhereRefIsNotBest(boot.sc3, refPredIndex=1, invertBayes=F) %>%
  as_tibble(rownames = "submission") %>%
  rename(bayes = value) 


# bayesFactorMatrix <- function(bootstrapMatrix) {
#   names <- colnames(bootstrapMatrix)
#   results <- purrr::map(1:ncol(bootstrapMatrix), function(ind) {
#     computeBayesFactorWhereRefIsNotBest(bootstrapMatrix, refPredIndex=ind, invertBayes=F) %>%
#       as.data.frame()
# }) %>% bind_cols
#   colnames(results) <- names
#   results
# }

---

Plot Results

Plot boxplot of all scores, coloring the boxes by Bayes factor.

plot_results <- function(results, bayes, subchallenge) {
  res <- results %>%
    as_tibble() %>%
    gather(submission, bs_score) %>%
    left_join(bayes) %>%
    mutate(bayes_category=case_when(
      bayes == 0 ~ "Reference",
      bayes<=3 ~ "<3",
      bayes>=3 & bayes <5 ~ "3-5",
      bayes>=5 & bayes <10 ~ "5-10",
      bayes>=10 ~ ">10")) 
  
  lvls <- colnames(results)
  res$submission <- factor(res$submission, levels = lvls)
  
  labs <- glue::glue("{res$submission}+")
  
  ggplot(res, aes(
      x=submission,
      y=bs_score,
      color=bayes_category
    )) +
    geom_boxplot() +
    theme_bw() +
    scale_color_manual(values = c(
      "Reference"="#d32e36", 
      '<3' = '#cf4d6f', 
      "3-5" = "#cc7e85",
      "5-10" = '#c5afa4', 
      ">10" = "#a8a6a4"),
      name = "Bayes Factor") +
    labs(x="Team", y=paste("Bootstrapped", subchallenge, "Score")) +
    theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))
  }

Subchallenge 1

SC1 ensemble results are a little odd. Adding in RYM has a pretty substantial effect on the mean ensemble.

plot_results(boot.sc1, bayes_top.sc1, "SC1")

NA
NA

ensemble_plot_data <-ensemble.sc1 %>% 
  pivot_longer(names_to = "model", values_to = "prediction", cols = -c("weight"))

lvls <- colnames(ensemble.sc1)[-2]
ensemble_plot_data$model <- factor(ensemble_plot_data$model, levels = lvls)
  
ggplot(data = ensemble_plot_data) +
  ggbeeswarm::geom_beeswarm(aes(x = model, y = prediction), priority = "density", size = 0.1, cex = 0.4) +
      theme(axis.text.x = element_text(angle = 90, hjust = 1, vjust = 0.5))

Subchallenge 2

SC2 results look more conventional : combine top models and increase performance a bit, but effect eventually disappears.

plot_results(boot.sc2, bayes_top.sc2, "SC2")

Subchallenge 3

SC3 results look more conventional : combine top models and increase performance a bit, but effect eventually disappears.

plot_results(boot.sc3, bayes_top.sc3, "SC3")